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Nonthermal fixed points of the dynamics of a dilute degenerate Bose gas far from thermal equilib- 
rium are analysed in two and three spatial dimensions. Universal power-law distributions, previously 
found within a nonperturbative quantum-field theoretical approach and recently shown to be re- 
lated to vortical dynamics and superfiuid turbulence [Phys. Rev. B 84, 020506(R) (2011)], are 
studied in detail. The results imply an interpretation of the scaling behavior in terms of indepen- 
dent vortex excitations of the superfiuid and show that the statistics of topological excitations can 
be described in the framework of wave turbulence. The particular scaling exponents observed in 
the single-particle momentum distributions are found to be consistent with irreversibility as well as 
conservation laws obeyed by the wave interactions. Moreover, long-wavelength acoustic excitations 
of the vortex-bearing condensate, driven by vortex annihilations, are found to follow a nonthermal 
power law. Considering vortex correlations in a statistical model, the long-time departure from the 
nonthermal fixed point is related to vortex-antivortex pairing. The studied nonthermal fixed points 
are accessible in cold-gas experiments. The results shed light on fundamental aspects of superfiuid 
turbulence and have strong potential implications for related phenomena, e.g., in early- universe 
inflation or quark-gluon plasma dynamics. 
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I. INTRODUCTION 

Turbulence is a generic phenomenon observed in the 
relaxation dynamics of many-body systems far from ther- 
mal equilibrium [1] . It comprises a quasi-stationary flow 
of energy within certain incrtial regimes in momentum 
space [2] . Correlation functions like the energy spectrum 
and higher momenta of the velocity distribution exhibit 
universality and scaling [3, 4]. This, and quasistation- 
arity are the key characteristics rendering turbulence a 
nonthermal fixed point of the system's dynamics. 

In quantum many-body systems, from the formation 
of Bose-Einstein condensates in ultracold gases to quark- 
gluon plasmas produced in heavy-ion collisions and re- 
heating after early-universe inflation, nonequilibrium dy- 
namics governs many interesting phenomena. In this con- 
text, turbulence during thermalization is being studied 
with increasing effort [5-17]. A nonthermal fixed point 
of the evolution of a many-body system has the poten- 
tial to strongly affect the equilibration process by forcing 
the evolution to critically slow down before the system 
can thermalize. New scaling laws were found by analyz- 
ing non-perturbative quantum field dynamic equations 
[6-8, 10]. Analogous predictions for a dilute ultracold 
Bose gas were given in [9] , proposing strong matter- wave 
turbulence in the regime of long-range excitations. 

Superfiuid turbulence, also referred to as quantum tur- 
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bulence (QT) has been the subject of extensive studies in 
the context of helium [18, 19] and dilute Bose gases [20- 
23]. In contrast to eddies in classical fluids, vorticity in 
a supcrfluid is quantized [24, 25], and the creation and 
annihilation processes of quantized vortices are distinctly 
different [18, 19]. The observation of a Kolmogorov 5/3- 
law [3, 4] in experiments with supcrfluid helium [26-28] 
received much attention [29-35]. In particular, the role 
of the normal-fluid as compared to the supcrfluid com- 
ponent in the turbulent flow is under debate [18, 19]. 

Supcrfluid turbulence plays an important role in the 
context of the kinetics of condensation and the develop- 
ment of long-range order in a dilute Bose gas. This, as 
well as turbulence in its acoustic excitations has been dis- 
cussed in Refs. [36-42]. A possible observation of QT in 
ultracold atomic gases presently poses an exciting task 
for experiments [43-45]. Here we emphasize that the 
experimental study of superfiuid turbulence and, more 
generally, of nonthermal flxed points in ultracold Bose 
gases has strong potential implications for many other 
areas of physics. Besides vortical excitations this also in- 
cludes other (quasi-)topological excitations such as soli- 
tary waves in one-dimensional gases. 

A satisfactory ab-initio mathematical description of 
both quantum and classical turbulence is inherently dif- 
ficult due to the strong correlations building up within 
the system. Analytical results are known, however, in 
regimes where kinetic theory applies: In a dilute, degen- 
erate Bose gas the normal-fiuid component can vary at 
the expense or gain of the superfiuid part. As a conse- 
quence, the gas is compressible and so-called weak wave- 
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turbulence phenomena can occur for which scaling laws 
can be derived by analyzing kinetic equations [46-48]. 

Generically, however, the description in terms of wave 
kinetic equations such as the Quantum Boltzmann equa- 
tion breaks down in the infrared (IR) regime of long wave- 
lengths. For a Bose gas in this regime, single-particle oc- 
cupation numbers grow large and the description in terms 
of, for example, elastic two-to-two collisions becomes un- 
reliable. In the infrared limit, so-called strong wave tur- 
bulence is expected to occur. Recent developments pre- 
sented in Refs. [6-10, 14] allow one to set up a unify- 
ing description of scaling, both in the ultraviolet (UV) 
Quantum Boltzmann kinetic regime and in the infrared 
limit. In the IR regime, new scaling laws were found for 
a relativistic scalar field by analyzing non-perturbative 
Kadanoff-Baym dynamic equations with respect to non- 
thermal stationary solutions [6]. 

In this article, we present the details of our studies 
of the relaxation of two and three dimensional dilute 
Bose gases through stages of superfluid turbulence and 
the approach of a nonthermal fixed point, by means of 
simulations in the classical-wave limit of the underlying 
quantum field theory. In Section II, we summarize the 
quantum-field theoretical predictions of Refs. [9, 46], in 
particular for the scaling exponents of the single-particle 
momentum distribution and compare these with the nu- 
merically determined scaling. While we find excellent 
agreement, our results provide an interpretation of the 
nonthermal fixed points proposed in Refs. [6, 9] for the 
case of an ultracold Bose gas: The appearance of nonper- 
turbative infrared scaling reflects the presence of statisti- 
cally independent vortices, as previously reported in brief 
in Rcf. [16], sec also Figs. 2-5. This phenomenon appears 
to be distinctly different from the weak wave turbulence 
we observe in the UV part of the spectrum. 

In Section III, we describe a model of independent 
as well as pair-correlated vortical excitations [24] which 
allows the more refined interpretation of the scaling 
behavior during the different stages of the evolution 
presented in Section IV. We show that the stationary 
scaling is maintained by the presence of energy (UV) 
and particle (IR) fluxes, further supporting the ana- 
lytic theory. In Section V, we study the power-law 
distributions of the compressible and incompressible 
contributions to the flow pattern. This adds to the 
clear understanding of the bimodal scaling laws in the 
overall momentum spectrum. IR power-law spectra 
of subdominant compressible excitations suggest the 
presence of acoustic turbulence [46] on the top of the 
vorticity-bcaring quasicondcnsate. In comparison with 
recent experiments and analytical predictions, the 
velocity field probability distribution as well as the 
velocity statistics of individual vortices are studied. 
This observable is of great interest, since it has recently 
been used to experimentally verify a distinction between 
classical and quantum turbulence [49]. We finally show 
that the complete decay of the turbulent scaling is 
anticipated by the appearance of weaker IR power laws 



reflecting vortcx-antivortex pairing correlations. 



II. SUPERFLUID TURBULENCE AS A 
NONTHERMAL FIXED POINT 

We begin with a brief summary of the analytical and 
numerical results on dynamical fixed points and matter- 
wave turbulence reported in Refs. [9, 16] before we dis- 
cuss these results in terms of a statistical model of vortex 
excitations. Motivated by the original work presented 
in Rcf. [6] nonthermal fixed points of Kadanoff-Baym 
dynamic equations for time-dependent Green functions 
were analysed, for the case of an ultracold Bose gas, in 
Ref. [9]. At these fixed points of the dynamical evolu- 
tion of the system, the single-particle momentum distri- 
bution n(k, t) becomes (quasi-) stationary, n(k, t) — 0. It 
furthermore exhibits a characteristic universal power-law 
behavior, i.e., it scales according to 

n(sk) = s"«n(k) , (1) 

in a certain regime of momenta k. Here s is some posi- 
tive, real number and C a universal exponent which was 
determined from the dynamic equations for the field cor- 
relation functions. Different exponents resulted in dif- 
ferent momentum regimes. Exponents valid in the ul- 
traviolet regime of large ]k] were found to correspond 
to well-known fixed points of weak wave turbulence [46] 
in the respective systems. In addition to these, new, 
larger exponents were predicted in the infrared regime of 
small ]k] on the basis of a non-perturbative analysis of 
the Kadanoff-Baym dynamic equations derived from the 
two-particle-irreducible (2PI) effective action [6-10]. 

Simulations of the classical field equations for a rela- 
tivistic 0(Af)-symmetric scalar model, for A'' = 4, con- 
firmed the existence of the analytically derived scalings 
in the infrared regime [6]. As reported in Ref. [16], cor- 
responding simulations were performed for an ultracold, 
degenerate Bose gas in two and three spatial dimensions 
which demonstrated that for these systems the infrared 
scalings predicted in [9] rcfiect and are caused by the 
presence of quantized vortical excitations of the supcr- 
fiuid. In this way, the nonthermal fixed points can be 
related in a clear manner to topological excitations of 
the interacting coherent matter-wave field which shows 
that both, weak wave turbulence and macroscopic, topo- 
logical excitations of the field can be described within 
a unified ficld-thcorctic approach. Differently expressed, 
both, weak turbulent fiow and non-linear solitary bulk 
excitations are described in a unified manner as repre- 
senting a nonthermal critical fixed point of the system. 
In turn, the approach also implies that superfluid tur- 
bulence can be studied in a new way, in the frame of 
a universal quantum field theoretical approach. We re- 
mark that the relation between the new infrared expo- 
nents and topological excitations of the background (i.e., 
condensate) field gains further support by the topological 
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pattern formation described in Ref. [17] for a relativistic 
0(2)-synimetric scalar model, which serves as a partic- 
ular model for the reheating period after early-universe 
inflation. 

We will discuss in detail vortical excitations in a degen- 
erate Bose gas and their relation to the above nonthermal 
fixed points. Before this we summarize briefly the rele- 
vant results of Refs. [9, 16]. To be specific, throughout 
this article we consider an ultracold Bose gas of atoms 
in d = 2, 3 space dimensions interacting through s-wavc 
collisions which is described by the Hamiltonian 



H 
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where the time and space dependent fields $ = <i>(x, t) 
satisfy Bose commutation relations, and where the cou- 
pling g = Aira /m in three dimensions is defined in terms 
of the s-wave scattering length a. Here and in the fol- 
lowing we set h = I, see Appendix B for more details. 



A. Weak wave turbulence 

Suppose the generic case that for sufficiently large 
momenta |k| = k occupation numbers n(k, t) = 
(<I>^(k, t)$(k, t)) are small enough such that for a given 
coupling g the Quantum Boltzmann Equation (QBE) 

dtnk^I{k,t), (3) 

/(k,t)= J d'^d'^gd'V|rkpqr|''5(k + p-q-r) 

X [(^k + l)(np + l)nqnr 

- nkrip(riq + l)(nr + 1)], (4) 

describes the evolution of rik = n(k, t) under the effects 
of collisions. In our case, the transition matrix element 
squared |TpkqrP is a numerical constant proportional to 
g^ and thus independent of momenta. Zeroes of the scat- 
tering integral /(k) correspond to fixed points of the time 
evolution within the regime of applicability of the QBE 
[46] . Most prominent amongst these are the thermal fixed 
point corresponding to the system in thermal equilibrium 
and the trivial fixed point where the occupation number 
is independent of k. At both fixed points the scattering 
integral vanishes and n(k, t) becomes independent of t. 
Note that both, the trivial distribution and the thermal 
Bose-Einstein distribution in the Rayleigh- Jeans regime, 
for a;(k) ^ fc^, show a power- law behavior of the form (1) 
with C = and C = 2, respectively. 

The theory of weak wave turbulence [46] allows to ana- 
lytically derive further, nonthermal fixed points at which 
the occupation number ?i(k) obeys a scaling law of the 
form (1) and, in general, C ^ 2. As in classical turbu- 
lence of an incompressible fluid one assumes that univer- 
sal scaling appears within a certain regime of momenta, 
the inertial range. According to this picture, outside the 



scaling regime kinetic energy enters the system from an 
external source and/or is dissipated into heat, whereas 
there are no sources and sinks within the inertial inter- 
val. Instead, kinetic energy is transported from momen- 
tum shell to momentum shell without loss or gain. To a 
good approximation, this process is described by a conti- 
nuity equation in momentum space, with a momentum- 
independent, radially oriented current vector.^ A cen- 
tral aspect of weak-wave-turbulence theory is that the 
QBE can be cast into different such equations [46], for 
the radial number density N{k) = (2fc)'^~^7rn(A;) and the 
energy density E{k) = (2fc)''-i7r£(fc), e{k) = u}{k)n{k), 



dtN{k,t)^-dkQ{k), 
dtE{k,t) ^ -dkP{k). 



(5) 
(6) 



Depending on whether the radial particle current Q{k) ~ 
{2kY-^iTQk{k) or energy current P{k) = {2k)'^-^nPkik) 
is taken to be independent of k, one derives different 
scaling exponents. The resulting exponents are 
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These exponents can be obtained by simple power count- 
ing: Combining Eqs. (3) and (5) gives the radial relation 
dkQ{k) ~ k'^~^I{k) which implies that stationarity re- 
quires k'^I{k) to become fc-indcpendcnt, i.e. scale like 
Counting all powers of k in I{k), Eq. (4), in the wave- 
kinetic regime where the terms of third order in the oc- 
cupation numbers dominate the scattering integral, this 
requires n{k) ~ Analogously one infers the ex- 

ponent (p^ from the balance equation (6) for the energy 
density e(fc) ~ k'^n{k). Despite this simple procedure, 
the existence of the respective scaling solutions has to 
and can be derived rigorously from the QBE by means 
of Zakharov conformal integral transforms [46] . We note 
that, similar to classical turbulence, the case d = 2 is 
special, where the scaling exponent Cp^ equals that of a 
thermal distribution in the Rayleigh- Jeans regime. 



B. Strong wave turbulence 

Given a positive scaling exponent C momentum occu- 
pation numbers nik) ^ grow large in the IR regime 
of small k. As a consequence, for a given coupling g, the 
QBE fails in this regime, where contributions to the scat- 
tering integral lik) which are of higher order than g'^n^ 
become important. 

To find scaling solutions in the IR, an approach be- 
yond kinetic theory is required. This is available through 
quantum-field dynamic equations derived from the two- 
particle irreducible (2PI) effective action or ^-functional 



^ Note that justification of tfiis assumption, i.e., locality of the 
transport, needs to be checked for each particular wavc-turbulent 
solution [46]. 
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ergy flow (P) in d dimensions were predicted in Ref. [9] 
to be 



(c) 



FIG. 1: (Color online) 2PI diagrams of the loop expansion of 
T2 [G] . (a) The two lowest-order diagrams of the loop expan- 
sion which lead to the Quantum-Boltzmann equation. Black 
dots represent the bare vertex ~ gS{x — y), solid lines the 
propagator G{x,y). (b) Diagram representing the resumma- 
tion approximation which, in the IR, replaces the diagrams 
in (a) and gives rise to the scaling of the T-matrix in the IR 
regime, (c) The wiggly line is the scalar propagator which is 
represented as a sum of bubble-chain diagrams. See text for 
more details. 



[50-52] beyond the 2-loop order in the expansion of the 
self-energy. See Ref. [9] for details of the procedure sum- 
marized in the following. The 2PI equations include 
the Dyson equation for the time-ordered Green function 
G{x,y) = {T^^x)^{y)) (here, we use four-vector nota- 
tion X = (xojx)), from which a time evolution equation 
(3) for n(k) is derived. As before one considers zeros of 
the scattering integral which in the dynamic theory reads 



(8) 



Here, k = (w, k), and p and F are the spectral and sta- 
tistical components of G, respectively, defined in coor- 
dinate space by F{x,y) = {{^'<{x),<i>{y)})/2, p{x,y) = 
i{[<^Hx),^y)]), G{x,y) = F{x,y) ~ {i/2)sgn{xo - 
yo)p{x,y). The corresponding contributions to the self 
energy 'E,{x,y) ~ 2iST2/SG{x,y) are defined in terms of 
G through a loop expansion of the 2PI effective action, 
sec Fig. 1. Rcsumming an infinite set of such loop di- 
agrams contributing to the 2PI effective action [53, 55] 
leads to a non-perturbative, effectively renormalized cou- 
pling in the dynamic equations [6, 9]. 

To derive the scaling solutions of the dynamic equa- 
tions for Green functions G(w,k) one assumes sepa- 
rate scaling of its spectral and statistical components 
according to p{s'^uj,s'k) = s^^+''/9(aj, k), F{s^Lu,s'k) — 
s~'^~'^F{uj,'k), s > 0. Here, z is the dynamical scaling 
exponent accounting for a different scaling in lo as com- 
pared to k. The scaling exponent ^ is related to n by 
n(sk) = s^^^^''n(k). k is derived from the condition 
that the scattering integral (8) vanishes, making use of 
Zakharov transformations, while the anomalous scaling 
exponent rj remains undetermined by this and, for the 
first, can be assumed to vanish for the d = 2 and 3 di- 
mensional cases considered in this paper. The IR scaling 
exponents for radial quasiparticle fiow (Q) and radial en- 
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IR 



2, Cp = + 2 + ^ . 



(9) 



In situations where a quasiparticle picture applies the 
dynamical exponent z corresponds to the homogeneity 
index of lo: aj(sk) — s^aj(k). In this case the non- 
perturbatively rcsummcd effective coupling can be re- 
lated to the diagonal elements of an effective many- 
body T-matrix, T^^^^ = T£+p^q+r in the kinetic Boffz- 
mann formulation. In the scaling regimes this scales as 
\Tf\ ^ \T^%\ - \9Ck^-V[l + G'gk^-W]\, k = |k|, 
where G' is some constant which fine-tunes the position 
of the transition from UV to IR scaling: For small Uk 
and 2 = 2 one recovers the UV case discussed in the pre- 
vious section, i.e., T^^ is a constant independent of k. 
For large n^, the second term in the denominator dom- 
inates which, assuming scaling of rik ~ ^ ''7 implies a 
power-law behavior |T^^|^ ^ kHC-d+z) g^j^^^^ -^^ turn, the 
modified scaling (9) of in the infrared regime of small 
wave numbers as compared to the UV regime discussed 
before. Physically, the renormalized T-matrix implies a 
reduction of the effective interaction strength in the IR 
regime of strongly occupied modes [6] . As a consequence, 
single-particle occupation numbers rise, towards smaller 
wave numbers, in an even steeper way than in the weak- 
turbulence regime. 



C. Wave-turbulent scaling and vortices 
in a Bose gas 

We now briefly review the relation between wave- 
turbulent scaling and the appearance of vortical excita- 
tions in an ultracold degenerate Bose gas. In [16] the 
results of semiclassical simulations of the gas dynam- 
ics were reported, obtained by solving the classical field 
equation 



i9t(^(x, t) 



-|^+5l0(x,t)|^ 



(x,t) (10) 



in a box with periodic boundary conditions. See Ap- 
pendix B for details on the simulations and on lattice 
units. The initial field (/)(x, 0) was prepared by macro- 
scopically populating a few of the lowest momentum 
modes in the computation such that the resulting conden- 
sate density in configuration space varied between zero 
and some maximum value. Quantum noise is taken into 
account by adding a small random contribution to each 
field mode. Vortical excitations were created in large 
numbers, within shock-waves forming during the non- 
linear evolution of the coherent matter-wave field. Fig. 2 
shows an example of the condensate phase distribution 
over a numerical lattice in d = 2 dimensions, with vor- 
tices and antivorticcs indicated by white crosses and cir- 
cles, respectively. Fig. 4 depicts the location of vortex- 
ring cores in a d = 3-dimcnsional system. For videos 
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2r)(i 

Lattice site x 

FIG. 2: (Color online) The plot shows the phase of the com- 
plex field for a single run of the evolution in d = 2 dimensions. 
Parameters: = 3 x 10"^, N = 4x10^, Ns = 512, t = 46340, 
coordinates in lattice units of Ua, see App. B. White rings 
(crosses) mark vortices (antivortices). For a video of the evo- 
lution see [54]. 
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Radial momentum k 

FIG. 3: Single-particle mode occupation numbers are shown 
as functions of the radial momentum k at the time shown in 
Fig. 2. k is in lattice units, see App. B. Note the double- 
logarithmic scale. The development of a IR scaling with 
n{k) ^ fc~* coincides with the presence of superfluid vortices, 
with additional wave-turbulent or thermal background short- 
wavelength fluctuations in the UV. Dashed lines indicate the 
filling of the initially occupied modes with fc > 0. 



of the evolution see [54]. Both figures show a snapshot 
of the system at some time after the vortical excitations 
have formed. At these times, in d = 2, part of the vortex- 
antivortex pairs have re-annihilated vifith each other al- 
ready, and in c? = 3 some of the rings have undergone 
reconnections and also disappeared by shrinking to zero 
size. 

Runs were repeated many times (^ 100 times in d ~ 2 
and ~ 10 times in d = 3) for an ensemble of initial config- 
urations differing by statistical noise. The number of runs 
was chosen such that the statistical error arising from 




FIG. 4: Vortex line tangles are shown in three dimensions. 
Black dots indicate where the amplitude of the complex field 
falls below 5% of the mean density n. Parameters are: p = 
4 X 10"'', N = 6Ax 10^°, Ns = 512, t = 3276, coordinates in 
lattice units, see App. B. For a video of the evolution see [54]. 




0.03 0.1 0.3 

Radial momentum k 

FIG. 5: Single-particle mode occupation numbers are shown 
as functions of the radial momentum k for the snapshot in 
Fig. 4. k is in lattice units. Note the double-logarithmic scale. 
The development of a IR scaling n(fc) ~ k~^ coincides with 
the presence of superfiuid vortex lines. Dashed lines indicate 
the filling of the initially occupied modes with fe > 0. 



run-to-run fluctuations was reduced to a value on the or- 
der of the size of the symbols used in the figures. In the 
regime of momenta with large occupation numbers where 
quantum statistical fluctuations play little role, correla- 
tion functions like the spectrum could be computed 
in a quasi exact way by averaging at a given time over 
the ensemble. In this way, the fiow was analysed in terms 
of the ensemble- and angle-averaged single-particle spec- 
trum 



n{k) 



id-l 



(r(k)</'(k)) 



ensemble 5 



(11) 



as a function of radial momentum k = jkj. Figs. 3 and 5 
show this spectrum for d = 2 and 3, at the times chosen 
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in Figs. 2 and 4, respectively. Most importantly, these 
spectra show scaling, Eq. (1), within a range of momenta 
between about the maximum of the initially occupied mo- 
mentum modes indicated by the vertical lines and k ~ 0.3 
in lattice units. The exponent C is in agreement with the 
field theoretical prediction = d + 2 given in Eq. (7). 

It was found that the annihilation processes which 
eventually destroy the vortical structure while coherence 
builds up in the system are very slow and thus stabilize 
the scaling solution over a long time. In particular, only 
after the last vortcx-antivortex pair has annihilated and 
the last ring shrunk to zero the power law n[k) ~ A: ^<5 
breaks down and a thermal distribution of particles is 
left over. Characteristic times are t ~ 0(10^) for vortex 
formation, t ~ C(IO^) for stabilization of the scaling so- 
lution and t - 0(10") - 0(10'^) for the last vortices to 
annihilate. 




FIG. 6: (Color online) Left panel: Sketch of a random vortex- 
antivortex distribution underlying the IR scaling close to the 
nonthermal fixed point. Right panel: Correlated vortex dis- 
tribution causing a modification to weaker pair scaling in the 
IR, for momenta smaller than the inverse of the average pair- 
ing length. 

A. Independent vortices in d = 2 



III. POINT- AND LINE- VORTEX MODELS 

As we have seen above, the scaling of the momen- 
tum occupation numbers in the IR regime of large wave 
numbers, predicted within a non-perturbative analysis 
of strong wave-turbulence, correspond, for 2- and 3- 
dimensional degenerate Bose gases, to the appearance 
of macroscopic vortical excitations. In the following, we 
analyse the observed scaling spectra by comparing them 
to the single-particle momentum distributions for a set 
of randomly positioned point vortices (vortex lines) in 
d = 2 {d = 3) dimensions. We will show that uncor- 
related vortices (vortex lines) are sufficient to yield the 
infrared scaling with exponent C}^, see Eq. (9). We will 
also show that beyond this, pair correlations between vor- 
tices and antivortices as well as configurations with small 
rings well separated from each other can give rise to a 
further scaling exponent deviating from (see Fig. 6). 

The point (line) vortex model employed here was in- 
troduced by Onsagcr in 1949 [24]. It describes the com- 
plex flow pattern in terms of the statistical mechanics 
of interacting classical point objects. This model has 
been constructed discrete vorticity approximation 
of classical fluid turbulence, but it is even more suitable 
to describe supcrfluid turbulence consisting of quantized 
vortices. Further applications include plasma physics or 
stellar dynamics [56-58]. 

In the following, we study the spectrum of a set of vor- 
tices, flnding different scaling regimes dependent on vor- 
tex correlation functions. With this approach, nonther- 
mal fixed points in an ultracold Bose gas can be related 
to the statistics of vortices. Numerical calculations, sam- 
pling field configurations of static randomly positioned 
vortices, confirm the analytical result. Finally, the veloc- 
ity probability distribution of a vortcx-dominatcd fiow 
[49, 59] is derived and numerically confirmed. 



In two dimensions, an isolated, singly quantized vortex 
is described by the complex field (j)(r, ip) = ■\/?i(r)e*'^. As 
the r-dependence of the density n{r) only becomes im- 
portant at small scales on the order of the healing length 
^ — l/^/2mgn where our simulations are dominated by 
thermal excitations, we will omit it in the following and 
assume n to be uniform. 

A system of M vortices in d = 2 dimensions can be 
defined as (/)(x) = nf-'^0i(x), where (/)i(x) = (/)(x — x^) 
is the single vortex field centered around x^. We derive 
the corresponding particle spectrum by considering the 
hydrodynamic velocity field v = \Iip/m. Denoting the 
velocity field of a single vortex as v(x) (see Appendix 
Al) we can express the mean classical kinetic energy 
density of the velocity field as 

£;,(x) = |(|v(x)|2) 

= f (I / dVv(x-x')p(x') (12) 

where p(x) = X]f=i 'i^i<5(x — x^) defines the spatial distri- 
bution of vortices with winding number tii ~ ±1. Here 
and in the following, (•) denotes an ensemble average over 
different realizations of the classical field (?!)(x). We derive 
the low-A; scaling of n{k) from the kinetic-energy spec- 
trum E^(k), given by the angle- averaged Fourier trans- 
form of i?v(x), taking into account that at low fc, the 
single-particle spectrum is dominated by the superfluid 
velocity field v, i.e., 

n{k) ~ 2mk''^E^{k). (13) 

One has, from Eq. (12) 

£;.(!.) ~(|v(k)|2) = (|p(k)|2|v(k)|2), (14) 

with 

M 
i,3 



7 



Below the healing length scale = 2sin(7r/2^) (lattice 
units), the modulus of the velocity field of a single vor- 
tex scales as |v| ^ k^^ and is radially symmetric in mo- 
mentum space. Hence, the angle-averaged single-particle 
spectrum scales like 



n{k) = k~^ I M + 2'^K,KjJo{kkj) j . (16) 

Here, Jo{y) = (27r)~-^ /^^ d0cos(y cos((?)) denotes the 
zeroth-order Bessel function and Uj = |xi — Xj| is the 
distance between vortices i and j. 

Assuming that the positions x.^ of the vortices are un- 
correlated one can take the average over relative positions 
lij within the area Vr ~ ttR^, 



— / dllJoikl) = 2 
Jo 



Ji{kR) 
kR ' 



(17) 



and, for fixed A:, the limit i? — > oo. Hence, the second 
term in brackets in Eq. (16) vanishes and one finally ob- 
tains the scaling [60] 

n{k)r^k-^. (18) 



B. Independent vortex-antivortex pairs 

In a two-dimensional superfluid containing vortices 
and antivortices an effectively attractive force between 
the two species can lead to pair correlations. We study a 
signature of this feature in the single-particle momentum 
spectrum by applying the point vortex model introduced 
above to the case of vortex-antivortex pairs. 

As a first step we calculate the velocity field vva for 
a vortex-antivortex pair with the vortex situated at xi 
and the antivortex at — xi. 



VVA v(x - xi) - v(x + Xl). 



(19) 



The squared velocity field far away from the center of 
the pair can be obtained via a dipole approximation 



|x| ^ |xi| which yields the scaling vva 



■. Hence, in 



Fourier space, the pair velocity field scales as |vvAp ~ fc*' 
for low momenta. In this regime, the vortex-antivortex 
pair can again be treated as a point-like object with mod- 
ified velocity scaling. To obtain the infrared scaling of a 
set of random vortex-antivortex pairs, we define a spatial 
pair distribution /3pair(x) = ^(^ ~ Xi), with x,; denot- 
ing the center of the i-th vortex-antivortex pair. Then, 
the analysis performed for random vortices above can be 
adopted. Therefore, in the case of independently dis- 
tributed pairs, the approach predicts the infrared scaling 
of the occupation number to be the same as for a single 
pair, i.e. n{k) ~ 



C. Pair correlated vortices in d = 2 

The considerations from Sect. HI B are only valid in 
the far IR (or equivalently for pair size going to zero). 
In our numerical simulations we found, however, that 
for large evolution times, a finite minimum distance 
emerged between vortices and antivortices, and also be- 
tween vortices with the same circulation, see Sect. VD. 
To take into account this observation and to analyse 
the full spectrum, we go back to writing the distribu- 
tion p(x) = p^(x) — p^{x) as the sum of distributions 
/9^(x) = ^ ^Y) of M vortices and p^(x) = 

J2iLi ~ ) of antivortices. Hence, 



J d'xd'x' e'*'("-"''C(x,x'). 



(20) 



with C(X, X') = ( p^p^, ) - ( pVpA ^ _ ^ pApV ^ + ( ^ 

This allows for a derivation of the kinetic-energy distribu- 
tion in terms of correlation functions of vortex positions. 

We model pairing by the density-density correlation 
functions 



M 



^^V(A)^V(A)^ = i^5(x-x') + Px,xS 



Vr 

V(A)„A(V)^ = 

VrVx 



(21) 



0(A- |x-x'|) + Px,x'. (22) 



where V\ = ttA^ is the area in which the theta function 
equals one. The contributions 



M (M - 1) 
Vr{Vr - ^a) 



6i(|x-x'|-A) (23) 



take into account that, besides the pairing, vortices and 
antivortices avoid each other in the dilute gas, keeping 
a minimum distance A. This is due to vortex- vortex re- 
pulsion and fast vortex-antivortex annihilation on small 
distances. The functions Px.x' cancel out in Eq. (20).^ 

Inserting Eqs. (21), (22) into Eq. (20), evaluating the 
Fourier transform and angular averaging gives 



\\pik)\' 



2M 1 



kX 



Ji(fcA) 



(24) 



The expansion of the integral for fc ^ 27r/A yields the 
leading-order result 



M{kX)^/A + 0{k^) 



(25) 



and, from Eqs. (13) and (14), the same occupation- 
number scaling as for independent pairs. 



nik) 



(26) 



^ If different avoidance scales A apply for vortices and antivortices, 
the terms do not cancel, but the remaining term does not alter 
the results for pair scaling derived here. 



8 



IP 



o 

O 





n(k) . - 








fc-' 


- 


>■ fc-« - 




,1 , iS 



Radial momentum k 

FIG. 7; Radial momentum distribution for a randomly dis- 
tributed set of M = 80 bound vortex-antivortex pairs, on a 
= 1024^ grid. Note the double-logarithmic scale. The vor- 
tices are positioned according to the probability distribution 
(27). At low momenta the power law is consistent with the 
existence of random vortex pairs nt ~ k~^. Above fcpair, the 
distribution exhibits the scaling of an ensemble of indepen- 
dent vortices, ~ while at momenta larger than the 
healing-length scale one observes the vortex-core scaling 



At momenta k 3> 27r/A the independent- vortex scal- 
ing n{k) ^ fc^^ is restored. We note that the in- 
frared result (26) can also be achieved by choosing 
more general pair correlations of the form (PxPx') ^ 
(p^p^,) = M6[2y«7r(Af,_ - Af„iJ]-i|x- xT-^^(A„,ax - 
|x - x'|)6'(|x - x'l - Amin) + -Px,x', with £ > and 
Amin < Aniax, which is closer to our numerical obser- 
vations and includes the case of a fixed "binding length" 
Amin — >■ Amax- As before, pair scaling k~'^ is found for 
fcAniax ^ 1, critical scaling fc~^ for fcAmax ^ 1, irrespec- 
tive of e. We supplement our discussion of scaling in the 
point vortex model vifith numerical data obtained by av- 
eraging over an ensemble of field configurations in d = 2 
dimensions. These configurations were constructed by 
multiplying uncorrelated single-vortex fields centered at 
positions according to the probability distribution 



M 



■1 I •■ 



l[p('H.T)p^'H^Y (27) 



with pW(xr) = i/Vn,p'-^H^Y,^t) = v^'e{x-\^Y ~ 

xf I ) , where we neglect that unpaired vortices avoid each 
other, i.e., choose A = 0. They do not represent stable 
solutions of the classical field equation nor do they con- 
tain sound-wave and related excitations yet which would 
build up through the interactions of the field and the 
vortices. 

The resulting momentum spectrum is presented in 
Fig. 7. Three scaling regimes can be observed. At 
low momenta the power law is consistent with the ex- 
istence of random vortex pairs Uk '^-^ Above fcpair = 
2sin(7r/2A) ~ tt/A, choosing A = 25 (lattice units), the 



distribution exhibits the scaling of an ensemble of inde- 
pendent vortices, Uk ~ k~^, while at momenta larger 
than the healing-length scale k^ one observes the vortex- 
core scaling ^ k~^. The above results reflect that in a 
vortex dominated flow, particles with low momenta are 
found far away from vortex cores. In the case of pair- 
ing, the fiow field far away from the cores is given by the 
field of a vortex pair, and the low-momentum scaling fol- 
lows the pair-field scaling. Particles closer to the vortex 
cores pick up a higher momentum. Above fcpair the field 
will be dominated by the field of a single vortex. Note 
that Fig. 7 shows the result of a numerical calculation 
in which we have sampled field configurations of random 
static vortices. Dynamical simulations of the classical 
field equation Eq. (10) will be shown in Sect. IV. 

We close with recalling the picture Onsager developed 
in Ref. [24] of thermodynamic equilibrium states of a 
fixed number of vortices and antivortices in two dimen- 
sions. He used the Hamiltonian of vortical flow in two 
dimensions [61], 



1 



M 



ZTT ^ — ' 



(28) 



to describe the dynamics of a system of M vortices in a 
superfluid which hence behave like a Coulomb gas. Here, 
the position of the i-th vortex is denoted as = (cc^, y^). 
Due to the fact that the x and y coordinates of each 
vortex are canonical conjugates, phase space is identical 
with configuration space of the vortex positions. Hence, 
for vortices moving in a volume V the total phase space 
IS given by y*^. The Hamiltonian (28) implies that low- 
energy configurations feature vortices of opposite sign 
close to each other, whereas high-energy configurations 
require vortices of equal sign to group. Due to these con- 
straints, the number of configurations D,{E) available for 
the system at a given energy E decreases towards high 
and low energies, with a maximum at some intermediate 
E = Eq. According to Boltzmann, the entropy is S{E) = 
kBln{^{E)) and the inverse temperature 1/T ~ dS/dE 
is positive for E < Eq and negative for E > Eq. It fol- 
lows that positive-temperature states are characterized 
by vortex-antivortex pairing, while negative-temperature 
states feature vortices of the same circulation to clus- 
ter. At the point of maximum entropy S{Eq) and in- 
finite temperature, Onsager expected a state of uncor- 
related vortices and antivortices. This understanding of 
the nonthermal fixed point as a quasi-equilibrium state 
of vortices corroborates our result that the IR scaling be- 
havior at the fixed point corresponds to the appearance 
of topological excitations. 



D. Vortex loops in d = 3 

We now consider the three-dimensional case of vortex 
lines and loops. A formulation similar to the Onsager 
point vortex model in Sect. HI A is possible [34, 62, 63]. 
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FIG. 8: Radial momentum distribution of an elliptical vortex 
loop on a = 1024"^ grid. Note the double- logarithmic 
scale. The major and minor radius scales ka, for ra = 40, 
and kb, for rf, = 4, respectively, are indicated. These scales 
separate pair scaling n(fc) ~ k~^ as for a near-circular vortex 
ring, scaling n{k) ^ k~^ for two anti-circulating vortex lines, 
and n(fc) ~ fe""", as for a vortex ring or of a pair of straight 
vortex lines, corresponding to the scaling exponent ^g*^ at the 
nonthermal fixed point. See Sect. HID for more details. 



See also Refs. [64-67] for an extension including the dy- 
namics of tangles of vortex lines. 

We write the classical kinetic energy spectrum of the 
velocity field in terms of the vorticity density 



u;(x) = mV X v{x) 



Sv(k) 



(l^(k)l^ 



(29) 



(30) 



The vorticity is vanishing everywhere but on the vortex 
lines, i.e., assuming M individual vortex loops. 



drsKr)5(s.(r)-x). 

- Jo 



(31) 



In this expression, the vortex filaments are represented 
by the connected curves Si(T), parametrized by the one- 
dimensional coordinate r G {0,Li}, Li being the arc 
length of filament i. The loops are closed, s,;(Li) = 8^(0), 
possibly also accross the walls of the 3-dimensional vol- 
ume in accordance with periodic boundary conditions. 
s^(t) is the tangent vector along the filament at Si(r), of 
unit length |s^(t)| = 1. We parametrize the i-th vortex 
loop in terms of a single center coordinate and a relative 
curve, Si(T) = ~f Viir) and write the vorticity as 

'^W=E j d^y%-R.) drr^(r) 

x^(r,(r) -x + y). (32) 



M 



;(k) = 5]e^'^^'c:>,(k). 



(33) 



where ti)i(k) is the Fourier transform of the vorticity of 
the z-th vortex loop, given by 



^.(x)= / drrUT)J(r,(r)-x). (34) 
Jo 

The ensemble-averaged vorticity (|a;(k)p) becomes 

M 

(|u;(k)p) = ^(e^''(^'-^^)^,(k)c:>,(k)) . (35) 



Assuming that the shapes of the individual vor- 
tex loops are statistically independent of their po- 
sition, it follows that (e*'^(^'-^j'a;j(k)a>j (k)) = 
(e'*'^^'~^^)) (w,(k)a;j (k)). If the loops are also un- 
correlated among themselves, then (a;i(k)a;j(k)) ~ 
(|a>i(k)p)5ij. Hence, for statistically identical loops. 



(|u;(k)p)=A./(|cI,(k)p), 



(36) 



which means that the vorticity spectrum scales in the 
same way as the average vorticity of a vortex loop cen- 
tered at the origin. Finally, the scaling of the momentum 
spectrum n{k) = 2mk^'^E^{k) follows from that of the 
angle-averaged vorticity. 



n{k) 



dl7fc(|a;(k)|^ 



(37) 



For the case of two straight parallel vortex lines of op- 
posite circulation Eq. (36) is evaluated in Appendix A 2, 
for an ensemble of such paired lines in Appendix A 3. 
Making use of the procedure developed there, the case 
of circular vortex rings of radius r is discussed in Ap- 
pendix A 4. For the latter, the resulting angle-averaged 
momentum spectrum scales like n{k) ^ for momenta 
k <^ kj. = 2sin(7r/2r') and n(fc) ~ k~^ for momenta 
k ^ k,.. To include effects from squeezed vortex loops, 
we consider elliptical filaments. The ellipse is defined by 
a major radius and minor radius r?,. In Fig. 8, the 
angle-averaged momentum spectrum of an elliptical vor- 
tex loop is shown. Three scaling regimes can be distin- 
guished. For the lowest momenta, one has n{k) 
which equals the infrared scaling for a vortex ring. For 
momenta ka k <^ fc;,, one finds n{k) ^ k~^, which 
coincides with the infrared scaling of two anti-circulating 
vortex lines (see App. A 2). For the eUipse, above kt, the 
momentum distribution scales like n{k) ^ k~^, which 
is the same as the high-momentum scaling of a vortex 
ring (see App. A 4) or of a pair of straight vortex lines 
(App. A3). 
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FIG. 9: Velocity distribution (39) as obtained from a random 
distribution of M — 80 bound vortex-antivortex pairs, on a 
= 1024^ grid, as in Fig. 7. Note the double-logarithmic 
scale. At low momenta the velocity field follows the scaling of 
vortex-antivortex pairs (a = 2), whereas at higher momenta 
the distribution reflects the presence of vortices (q = 1). 
Note, that there is a deviation from a = 2 in the regime 
of low velocities. From simulations with higher resolution, 
we find the best agreement for a — 1.8. We remark, that in 
the limit M — > oo the velocity distribution is expected to be 
Gaussian [681. 
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FIG. 10: Single-particle mode occupation numbers as func- 
tions of the radial momentum k, for four different times of 
a run in d = 2 dimensions. Parameters are: g = 3 x 10~^, 
iV = 4 X 10**, iVs = 512. Note the double- logarithmic scale. 
An early development of a scaling n{k) ~ k~^ is followed by 
a quasi-stationary period of bimodal scaling with n{k) ~ fc~* 
in the IR, due to the presence of vortices, and n{k) ~ in 
the UV, corresponding to weak wave turbulence or thermal 
equilibrium. 



E. Velocity distribution 

The vortex model can provide insight into another ob- 
servable accessible in our numerical simulations, the ve- 
locity distribution. See Refs. [59, 68-70] for discussions 
of the velocity distribution and Ref. [49] for recent ex- 
perimental results. As discussed above, the velocity field 
scales, far away from any core, as ^ r~", with a = 1 for 
a single vortex and a = 2 for a vortex-antivortex pair 
in two dimensions. The velocity probability distribution 
7'(v) is calculated as 'P(v) = \dx/dv\ p{x), with spatial 
vortex distribution function p(x). For a uniform distri- 
bution p(x) = const., it follows after angular averaging 
that 

p(v) = |v|"2("+i)/". (38) 

It is numerically convenient to calculate the probability 
density of a single component of the field, e.g. 

= I d^-,7'(v)^^;^2("+l)/^ (39) 

In Fig. 9, the corresponding velocity distribution V{vx) is 
shown. In accordance with the analytical predictions two 
scaling regimes are observed. At low momenta the ve- 
locity field follows the scaling of vortex-antivortex pairs, 
whereas at higher momenta the distribution reflects a 
random distribution of vortices. 



IV. BOSE GAS APPROACHING THE 
NONTHERMAL FIXED POINT 

In this section we return to the process of the for- 
mation of vortical excitations and of the Bose gas ap- 
proaching the nonthermal fixed point characterized by 
the wave-turbulent scaling solutions discussed in Sect. II 
and Rcfs. [9, 16]. Thereby we first focus on the evolution 
of the single-particle momentum distributions and the 
emergence of nonthermal power laws. We relate these to 
structure formation in the form of vortical excitations. In 
order to identify in a clearer way the relevant processes 
in approaching the fixed point we compute the fluxes in 
momentum space. Remarkably, we find that the appear- 
ance of the particle scaling exponent C}q in the IR and of 
the energy exponent in the UV are compatible with 
predictions based on general arguments of irreversibility 
as well as energy and number conservation in a collision 
process [46]. 



A. Time evolution of the single-particle spectrum 

As discussed in Sect. II C, vortical excitations can be 
created in large numbers, within shock waves forming 
during the non-linear evolution of the coherent matter- 
wave field. Figs. 10 and 11 show the angle- and ensemble- 
averaged radial momentum spectra, Eq. (11), for a Bose 
gas in a box with periodic boundary conditions, in d = 2 
and d = 3 dimensions, respectively. Four snapshots are 
shown, taken at the dimcnsionless times t as indicated 
in each panel. See Appendix B for details on the sim- 
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FIG. 11: Single-particle mode occupation numbers as func- 
tions of the radial momentum k, for the four different times 
of a run in d = 3 dimensions. Parameters are: — 4 x 10"'*, 
iV = 8 X 10^, iVs = 256. Note the double-logarithmic scale. 
A period of bimodal scaling n{k) ^ (vortex lines, IR) 
and n(k) ~ k~'^ (weak wave turbulence, UV) is followed by 
trimodal scaling, which also exhibits pairing, i.e., n(fc) ~ 
in the far IR, induced by a set of far separated small vortex 
rings, in addition to n(fc) ~ k~^ (IR-): and thermal scaling 
n{k) ~ fc-^ in the UV. 



ulations and on lattice units. The initial field configu- 
rations were prepared by macroscopically populating a 
few of the lowest momentum modes in the computation 
such that the resulting condensate density in configura- 
tion space varied between zero and some maximum value. 
In Figs. 10 and 11, the panels representing the earliest 
times show the system after a brief initial evolution dur- 
ing which momentum is rapidly transported from the few 
initially occupied modes near zero over a comparatively 
large range of wave numbers. 

Shortly after this, vortical excitations are created 
which immediately causes the spectrum in d = 2 dimen- 
sions to exhibit a power-law behavior 2.85 ^ C ^ 3.0 
within a range of momenta k e [0.04 : 0.4], see the 
upper right panel of Figs. 10. Subsequently, the evo- 
lution slows down and a quasi-stationary period is en- 
tered. During an intermediate stage (bottom left panel 
of Fig. 10 and upper right panel of Fig. 11) of the vortex- 
bearing phase two distinct power laws develop which are 
in excellent agreement with the analytical prediction in 
Eqs. (7) and (9). While in the ultraviolet the expo- 
nent = d exhibits weak wave turbulence, Eq. (7), 
in the infrared, the exponent confirms the field theory 
prediction Cq* — d + 2, cf. Eq. (9). More specifically, 
in two dimensions, bX t — 5792, one observes scaling 
exponents 3.8 ^ C ^ 4.0 within a range of momenta 
k e [0.02 : 0.2] and 2.0 < C ^ 2.3 within a range 
of momenta k e [0.2 : 0.7], see the lower left panel 
of Figs. 10. At t = 262144, 4.0 < C < 4.2 within a 
range of momenta k E [0.02 : 0.2], see the lower right 
panel of Figs. 10. In d = 3 dimensions, at t = 820, 



one observes 4.8 ^ C ^ 5.0 within a range of momenta 
k e [0.08 : 0.4] and 3.0 < C < 3.1 within a range of 
momenta k G [0.5 : 1.7], see the upper right panel of 
Figs. 11. At ? = 1640, 5.0 < C < 5.1 within a range 
of momenta k € [0.05 : 0.5], see the lower left panel of 
Figs. 11. The appearance of the bimodal power laws cor- 
roborates results for a relativistic 0{N) model reported 
in Refs. [6, 10]. 

During the ensuing evolution, the weak-wave- 
turbulence scaling decays towards C = 2, reflecting a 
thermal UV tail. Note that in c? = 2, the weak-turbulence 
exponent C,p^ = 2 is identical to that in thermal equi- 
librium in the Rayleigh- Jeans regime, n{k) ^ T/k"^ 
[46]. In d = 3 we observe, at late times, a change of 
the infrared scaling behavior from C = (i+2 = 5to 

= 3, pointing to the development of pairing correla- 
tions, cf. Sects. IIIC, HID, and Sect. VD below. More 
specifically, dXt~ 26214, one observes a scaling exponent 
2-9 ^ C ^ 3.0 within a range of momenta k £ [0.03 : 0.1], 
see the lower right panel of Figs. 11. 

At late times, after the last vortical excitations have 
disappeared, we observe the entire spectrum to become 
thermal, i.e., exhibit Rayleigh- Jeans scaling with C = 2 
(not shown). We emphasize that thermal scaling of the 
single-particle occupation number with C = 2 applies de- 
spite the fact that quasiparticles with a linear dispersion 
are expected to thcrmalize in the regime of wave num- 
bers smaller than the inverse healing length 1/^ ~ 0.45 
{d = 2) and 1/^ ~ 0.22 {d 2). In the Bogoliubov 
approximation this is seen by taking into account the 
power-law dependence of the coefficients u| ^ ^ k~^ 
which contribute to n{k) ~ (u| + vl)T/k ^Tjk-"^ [9]. 

We are going to study more details of the time evo- 
lution of the system in a forthcoming publication and 
focus in the remainder of this paper on properties of the 
stationary scaling distributions. An important question 
in this context is, why the system selects the particular 
exponents Cp^ = d and C}q = d -|- 2 from the set of four 
possible exponents given in Eqs. (7) and (9). For this, 
the fluxes underlying the stationary but non-equilibrium 
distributions are relevant. 



B. Fluxes 

The timeline of distributions shown in Figs. 10 and 11 
suggests that the evolution of the gas involves a transport 
process from intermediate momenta around 0.05 . . . 0.2 
both towards lower and higher wave numbers, building 
up a bimodal power-law distribution. To describe the 
character of these transport processes we plot, in Fig. 12 
for d = 2 (upper set of four panels) and d = 3 (lower set), 
the radial particle and kinetic-energy flux distributions 
Qfe and Pfe , respectively. Note that the radial particle flux 
density is multiplied by gn to have the same units as 
the energy flux density P{k). These flux densities are de- 
fined through the balance equations (5) and (6), respec- 
tively, with kinetic energy density Sk — nkk'^/2m. The 
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FIG. 12: (Color online) Kinetic-energy and particle fluxes 
in d = 2 (upper set) and d = 3 (lower set), for the four 
different snapshots of Figs. 10 and 11, respectively. Note the 
logarithmic fc-axes. The appearance of the bimodal scaling 
coincides with a positive kinetic energy flux in the UV, and 
a negative particle flux in the IR. Flux units: [P] = [gnQ] — 
(Waf+'')-\ cf. App. B. 



Rcf. [46], one can show on the general grounds of the 
Boltzmann 77-theorem that a necessary condition for a 
non-equilibrium stationary distribution in our systems 
is energy damping in the region of large k. Moreover, 
energy and particle number conservation in the interac- 
tion of different momentum modes can be shovirn to im- 
ply the existence of at least one more sink, i.e., a region 
where the right-hand sides of Eq. (5) effectively has an 
additional damping term ^ T(k)n(k), with negative T. 
In between these sinks, a source region supplies the in- 
put to the bidirectional flux pattern towards the UV and 
IR. It was shown, moreover, in Ref. [46], within wave- 
kinetic theory, that generally, a positive, fc-independent 
flux transports energy, P > while a negative flux trans- 
fers particles, Q < 0. The only exception is the case of 
d — 2 where the thermal and the weak-wave-turbulence 
exponent can not be distinguished. Remarkably, this 
pattern remains valid in our case, beyond the UV weak- 
wave-turbulence regime, in the IR region where the ex- 
ponent emerges from a fixed point of the full dynamic 
equations. As already pointed out in Ref. [9], however, 
the derivation of the IR exponent Cq* ^ d+2 requires the 
existence of a sufficiently well defined quasi-particle dis- 
persion relation, suggesting a treatment in terms of the 
Quantum Boltzmann equation with a momentum depen- 
dent scattering matrix element to be applicable. 

From this point of view, the negative flux Q and scaling 
in the IR and the positive flux P and weak wave turbu- 
lence in the UV, as observed in the numerics, emerge as a 
necessary consequence of conservation laws and transport 
processes described by wave-kinetic transport equations. 



V. VORTICES, ACOUSTIC TURBULENCE AND 
THE DEPARTURE FROM THE FIXED POINT 



initial stage is governed by an infrared particle transport 
towards larger fc, also causing a flux gnQ of interaction 
energy. At intermediate times, Q{k) changes sign. This 
is accompanied by a positive kinetic energy transport in 
the UV, as observed in two-dimensional simulations in 
Rcf. [71]. 

We emphasize that the negative particle flux in the 
IR and the positive kinetic-energy flux in the UV co- 
incide with the appearance of the bimodal momentum 
distributions in Figs. 10 and 11 (bottom left panels). Al- 
though the derivation of the IR exponents requires the 
full dynamical theory with non-perturbatively resummed 
self-energies, the signs of the fluxes correspond to the re- 
spective scaling exponents, i.e., Cq^ in the IR, and (p^ 
in the UV, cf. Sect. II A. Moreover, at late times, the 
kinetic-energy flux P vanishes due to a thcrmalizcd UV 
momentum distribution, but Q still reshuffles particles 
and therefore interaction energy, keeping the system out 
of equilibrium close to the nonthermal flxed point. 

Wc finally remark that, as discussed in detail in 



In this last section we analyse the single-particle mo- 
mentum spectra obtained in our numerical simulations 
in view of their interpretation in terms of vortical ex- 
citations and wave-turbulence as implied by the point- 
and line-vortex models introduced in Sect. III. For this 
we first decompose the flow pattern of the system which 
has closely approached the nonthermal fixed point, into 
transverse (incompressible) and longitudinal (compress- 
ible) contributions. In this way we can show that the IR 
scaling is dominated by the incompressible part while in 
the UV the particles mainly belong to the compressible 
as well as a quantum pressure components. We find a 
further scaling exponent C ~ d -I- 1 for the subdominant 
compressible component in the IR which is interpreted 
as a signature of acoustic turbulence. 

Moreover, going beyond the results presented in [16] 
we identify signs of pair formation in the final stage. 
The scaling caused by this goes beyond the predictions 
from dynamical quantum field theory summarized in Sec- 
tion II, and we interpret it as a signature for the system 
leaving the fixed point again for final thermalization. 
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FIG. 13: (Color online) Occupation numbers as functions of 
the radial momentum k, as defined in Eqs. (40)-(43): Total 
single-particle occupation number n{k) (black dots), incom- 
pressible (solenoidal-flow) component ni(fc) (red circles) , com- 
pressible (rotationless) component nc{k) (blue filled squares), 
quantum-pressure component nq(fc) (grey open squares). Pa- 
rameters are the same as in Fig. 10 (lower right panel), for 
the run in d = 2 dimensions at the time t — 262144. Note the 
double-logarithmic scale. A scaling with /c~* ®® corresponds 
to a power-law exponent 5/3 for the kinetic energy in d = 2. 
See text for more details on the scaling exponents. Inset: 
Phase angle (^(x, t) as in Fig. 2. 

A. Superfluid turbulence and statistics of vortices 

To exhibit vortical flow and define the decom- 
position we use the polar representation (/)(x, t) = 
\/n{x, t) exp{i(/?(x, t)} of the field in terms of the den- 
sity n{'x.,t) and a phase angle (p{'x,t). This allows to 
express the particle current j ~ i{4>*\/4>~ 0V0*)/2 = nv 
in terms of the velocity field v = \7(p. 

With this, we decompose the kinetic-energy spectrum 
following Ref. [29, 30], splitting the total kinetic energy 
Ekin = /d''x(|V(/)(x,t)|2)/(2m) as E^in = E^ + into 
a 'classical' part 

E^ ^ ^ J d'^x ilV^^f) (40) 

and a 'quantum-pressure' component 

E., = ^ I d'^x(lVV^p). (41) 
2m J 

The radial energy spectra for these fractions involve the 
Fourier transform of the generalised velocities Wv = v^'^ 
and Wq = V-y/n, 

Es{k)^^J k''~^A^a{\^&m''), '5 = v,q. (42) 

Following Ref. [29, 30], the velocity Wv, which due to 
the multiplication of v with the density n becomes reg- 
ularized and can be Fourier transformed, is furthermore 
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FIG. 14: (Color online) Occupation numbers as functions of 
the radial momentum k, as defined in Eqs. (40)-(43) (see 
caption of Fig. 13 for more information). Parameters are the 
same as in Fig. 11 (lower left panel), for the run in d = 3 di- 
mensions at the time t = 1640. Note the double- logarithmic 
scale. A scaling with fc"^"® corresponds to a power-law ex- 
ponent 5/3 for the kinetic energy in d = 3. Inset: The black 
spots mark points around the vortex line cores where the den- 
sity falls below 5% of the mean density n. 

decomposed into 'incompressible' (divergence free) and 
'compressible' (solenoidal) parts, Wv = Wi + Wc, with 
V • Wi = 0, V X Wc = 0, to distinguish vortical superfluid 
and rotationless motion of the fluid. For comparison of 
the kinetic-energy spectrum with the single-particle spec- 
tra n(fc), we determine occupation numbers correspond- 
ing to the different energy fractions as 

ns{k) = k-''-^Es{k), 5e{i,c,q}. (43) 

The resulting spectra ni{k), nc{k), and nq(fc) add up to 
ns{k) = ni(fc) -f nc{k) + nq{k), which agrees with the 
single-particle spectrum up to small corrections, see Ap- 
pendix C. 

In Figs. 13 and 14, we depict the momentum distribu- 
tions of the occupation numbers ni{k), nc(fc), and nq{k), 
together with the previously shown total single-particle 
spectrum n(fc), each at a late time when the system is 
close to the fixed point, i.e., shows the predicted scaling 
both in the IR and the UV. Red circles denote nj, filled 
blue squares ric, and open grey squares riq. Qualitatively, 
the results are similar for d = 2 and d = 3. 

In the range of large wavenumbers, the spectrum is 
dominated by the compressible and quantum-pressure 
components. The scaling of these excitations exhibits the 
weak-wavc-turbulence exponent (jP^ = d. For smaller 
wave numbers the scaling changes to n{k) ^ k^'^^"^. The 
decomposition into the various components now shows 
that this switching to a different scaling in the IR is 
clearly due to the take-over of a different character of 
the excitations with a modified flow pattern accounted 
for by the incompressible (solenoidal) component of Wy. 
The fact that in this regime contributions from vortical 
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flow rij dominate is in accordance with our interpreta- 
tion of the strong IR scahng by a model of an ensemble 
of vortical excitations. As we have discussed in the pre- 
vious section the analytically predicted infrared power 
laws n{k) ~ k~'^~^ are consistent with a finite density of 
independent vortices and antivortices {d = 2) or vortex 
lines (d = 3). 

Moreover, we find that the compressible component in 
the IR represents the second strongest contribution to 
the flow and develops a nonthermal scaling as ^ fc"'^"-'^. 
We will discuss this result in more detail in Sect. VB 
below. 

For momenta larger than about k ~ 0.2 for d = 2 and 
k ~ 0.4 for d = 3 the compressible and quantum-pressure 
components dominate the momentum distributions. In 
the regime of intermediate momenta, above the scale 
ki ~ 27r/Z with I the mean inter- vortex spacing, where 
the incompressible flow and the rest are of roughly equal 
strength one finds a scaling of approximately n{{k) 
^-d-i-.5/3^ corresponding to Ei{k) ^ k~^/'^ for both the 
two- and the three-dimensional case. Scaling with this 
exponent, which is the same as in classical Kolmogorov 
turbulence in an incompressible fluid, has been found in 
simulations of superfluid turbulence in d = 2,3 before 
[20, 21, 29-34, 72]. In some of these cases, turbulent 
flow appeared in the time evolution of a system starting 
from a Taylor-Green configuration of vortex line tangles. 
We emphasize that our configuration after the creation 
of vortical excitations resembles more a state of the kind 
usually termed chaotic turbulence. It is unclear whether 
chaotic turbulence can be related to classical Kolmogorov 
turbulence as this does not bear near-classical bundles of 
equal-orientation vortex lines. Nonetheless we observe 
such 5/3 scaling in the range of intermediate momenta 
where the interplay between the incompressible and com- 
pressible components is most pronounced. 



B. Acoustic turbulence 

The IR scaling ~ fc^''^^ of the compressible compo- 
nent (blue filled squares) in Figs. 13, 14, and also Fig. 18 
below suggests an interpretation in terms of acoustic tur- 
bulence [37, 38, 42, 46] and corroborates the numerical 
findings reported in Ref. [73]. The scaling is persistent 
until late times, see, e.g.. Fig. 18, but decays after the 
vortices have disappeared. To check that the nonthermal 
scaling of the incompressible component is not an artefact 
of the decomposition of the vortical flow into compress- 
ible and incompressible parts, we show the momentum 
spectrum for an average over selected runs in d = 3, see 
Fig. 15 (lower panel), where the snapshot is taken shortly 
after the last vortex ring has disappeared. One can ob- 
serve that acoustic turbulence can survive for a limited 
period. At = 0(100), beyond the time when vortical ex- 
citations are negligible. In the case of d = 2, the effect 
is weaker but still present, see Fig. 15 (upper panel), 
where we show the averaged spectrum shortly after the 
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FIG. 15: (Color online) Late-stage acoustic turbulence: Occu- 
pation numbers defined in Eqs. (40)-(43) are shown as func- 
tions of the radial momentum k, from an average over ~ 10 
single runs in d = 2 dimensions (upper panel) and d = 3 
(lower panel), at times t = 0{lQ^) (d = 2) and t = 0(10'') 
(d = 3), triggered on the decay of the last vortical excita- 
tion. See caption of Fig. 13 for more information. Param- 
eters are: d = 2: 5 = 3 x 10"'^ , = 4 x 10®, Ns = 512; 
d = 3: g = 4 X 10"*, Af = 10^ iVs = 128. Note the 
double-logarithmic scale. The figure shows that, shortly after 
the last vortex ring has disappeared (incompressible compo- 
nent breaks down), compressible excitations exhibiting acous- 
tic turbulence scaling ~ fc"'''^^' remain present. 

last vortex-antivortex annihilation. Our interpretation of 
this flnding is that acoustic turbulence coexists with the 
dominant vortical flow. It is driven by vortex motion and 
especially vortex annihilation processes, which are known 
to produce compressible excitations. In d = 2, signa- 
tures of acoustic turbulence are less pronounced, which 
we attribute to reduced vortex dynamics or, equivalently, 
weaker driving of compressible excitations. 

C. Vortex velocities 

We furthermore investigated the velocity-field proba- 
bility distribution in the turbulent regime to compare our 
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FIG. 16: Velocity field probability distribution. Parameters 
are the same as in Fig. 11, for the run in d = 3 dimensions, 
at the time t — 1640. Note the double-logarithmic scale. The 
black line is a Gaussian fit to the data which shows that in the 
low-velocity regime, the distribution is dominated by Gaus- 
sian fluctuations. A high-velocity scaling with v^^ reflects the 
presence of uncorrelated vortices, cf. Eq. (39). 
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FIG. 17: (Color online) Comparison of the vortex velocity 
distribution with the velocity-fleld probability distribution in 
d = 2 dimensions at t = 32768. At this time an average of 
86 vortices are present in the system. Parameters: g = 3 x 
10-'5, N = 16x10^, Ns = 1024. Note the double-logarithmic 
scale. A high-velocity scaling with reflects the presence of 
uncorrelated vortices, cf. Eq. (39). Gaussian fluctuations are 
suppressed in the distribution of velocities of the individual 
vortices, measured through the motion tracking of the vortex 
cores. 
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FIG. 18: (Color online) Pairing effects: Occupation numbers 
are shown as functions of the radial momentum k. Parameters 
are the same as in Fig. 11, for the run in d = 3 dimensions, 
at the time t — 26214 (lower right panel). Note the double- 
logarithmic scale. The scaling ~ k~^ in the far IR reflects 
pair correlations present in far-separated small vortex loops. 
(See caption of Fig. 13 for more information). 



low momentum region. In d = 2 dimensions, see Fig. 17, 
we compare the velocity- field distribution with the veloc- 
ity distribution of the positions of the individual vortices. 
Similar to the three-dimensional case the IR part of the 
velocity-field distribution is characterized by Gaussian 
statistics, whereas the UV regime shows single-vortex 
scaling [59] . The vortex velocity distribution is obtained 
from the statistical analysis of a vortex-tracking algo- 
rithm which is designed to detect regions of low den- 
sity accompanied by a winding number equal to one. 
The method is analogous to the experimental set-up of 
Ref. [49], where particles are trapped inside vortex cores 
to study vortex velocities in d = 3. In agreement with 
the experiment. Fig. 17 shows that Gaussian random mo- 
tion is suppressed for vortices, but ~ scaling clearly 
persists. In this way, quantum turbulence can be dis- 
tinguished from classical turbulence, where a continuous 
distribution of vorticity favors a Gaussian velocity distri- 
bution. 



D. Pairing and departure from the fixed point 



data with the expectations summarized in Sect. Ill E. 
In d = 3 dimensions, see Fig. 16, the results corrob- 
orate theoretical and experimental results reported in 
Rcfs. [49, 59, 74, 75]. The low-momentum distribution 
is characterized by Gaussian statistics, whereas the UV 
regime clearly shows scaling V{vx) ~ predicted by 
the model of randomly distributed, uncorrelated vortex 
lines, see Eq. (39). Note that even at late times, when 
vortex rings shrink in size, no signs of vortex pairing, 
'P{vx) ~ ''^x'^: can be seen in the velocity distribution. 
This is due to dominating Gaussian fluctuations in the 



In Fig. 18 we show the decomposition of the single- 
particle spectrum into the previously defined incompress- 
ible, compressible, and quantum-pressure components, 
see Eq. (43), for the 3-dimensional system at i = 26214, 
as in the lower right panel of Fig. 11. At this time, only 
a small density of vortex rings has survived, which de- 
cays under the influence of the noise field, cf. also Refs. 
[76, 77]. During the period when only a few small vor- 
tex rings remain one can observe a decrease of the in- 
frared scaling exponents of the occupation number and 
its incompressible part to the value C = 3. Since vortical 
excitations are still dominating the spectrum, we can in- 
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FIG. 19: (Color online) Pairing effects: Occupation numbers 
are shown as functions of the radial momentum A; in d = 2 
dimensions, at the time t — 262144: Occupation number ob- 
tained from averaging over all runs (black dots); Occupation 
number obtained from averaging over selected runs featuring 
tightly bound vortex antivortex pairs with maximum pair cor- 
relation length A < 5 (red circles). Parameters: <; = 3 x 10""", 
N = 10^ , Ns ^ 256. Note the double-logarithmic scale. A 
scaling with reflects the presence of correlated vortex- 
antivortex pairs in d = 2. 



terpret this observation in terms of the statistical point 
vortex model introduced in Sect. Ill, as vortex-antivortex 
pair correlations, see Sect. HID. This is consistent with 
the snapshots of individual runs showing small vortex el- 
lipses (see inset of Fig. 18). The scaling transition can 
be identified with the scale kh = 27r/2r;,, with estimated 
minor radius r?, ~ 15. 

Fig. 19 shows the average occupation number spec- 
trum att = 262144 in c? = 2 dimensions, as in the lower 
right panel of Fig. 10. In the average over all generated 
runs pairing effects can hardly be observed. However, 
for selected runs with small maximum pair correlation 
length A < 5, as introduced in Sect. IIIC, the pair scal- 
ing appears. We find that in the selected runs (~ 0.2% 
of total number of runs), at the chosen point of time 
one or at most a few vortex-antivortex pairs with pair- 
ing length smaller than the distance between pairs are 
present. This constitutes the final period of the evolu- 
tion, shortly before the last pair has disappeared through 
mutual annihilation and the system fully thermalizes. 
The generic configuration in d = 2 dimensions during 
the preceding time interval of critical slowing down close 
to the nonthermal fixed point is characterized by ran- 
domly positioned vortices and antivorticcs, with pairing 
correlations nevertheless present. Pairing is seen in the 
density-density correlation function between vortices and 
antivortices shown in Fig. 20, from which we read off a 
largest scale Amax — 100 . . . 150, cf. Sect. IIIC. This im- 
plies pair scaling below k ~ 0.02... 0.03, not visible in 
the few (black) full circles in Fig. 19. We observe that, 
in d = 2, scattering between the vortices of equal and op- 
posite circulation can lead to the increase of the pairing 
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FIG. 20: (Color online) Normalized pairing (density-density) 
correlations of vortices and antivortices in d = 2. Parame- 
ters as in Fig. 10, at the time t = 131072. Red crosses: Pair 
correlations of vortices and antivortices within the same cir- 
culation. Blue stars: Pair correlations between vortices and 
antivortices. The rise above 1 indicates an effective "binding" . 
Black crosses: Resulting correlation function C introduced in 
Eq. (20). 

length between correlated vortices and antivortices. This 
is in contrast to the situation for d 3 where the scat- 
tering between vortex rings rarely leads to an increase 
in the size of the rings. It is energetically advantageous 
for the rings to shrink in size, and thus, at a particular 
late time of the evolution the presence of separated small 
rings is generic and "pairing" is seen in the spectrum. 

In summary, our analysis gives a picture of the neces- 
sary conditions and the character of the nonthermal fixed 
point. It also shows that pairing effects are a first signa- 
ture for the system leaving this fixed point again for the 
final move to thermal equilibrium. 



VI. SUMMARY 

We have studied in detail superfluid turbulence in two- 
and three-dimensional dilute Bose gases by means of 
simulations in the classical-wave limit of the underly- 
ing quantum field theory. A focus is set on the iden- 
tification and characterization of stationary scaling so- 
lutions, in particular for single-particle momentum dis- 
tributions. Characteristic exponents C of n{k) corrob- 
orate the analytical predictions made in Ref. [9]. Our 
findings suggest that local field expectation values and 
short- to intermediate-range coherence, including topo- 
logical excitations, are at the basis of the infrared power 
laws predicted within a full nonperturbativc dynamical 
field theory [6-10, 17] employing two-particle-irreducible 
effective action techniques. 

We have shown that the stationary scaling is main- 
tained by the presence of particle (IR) and energy (UV) 
fluxes, originating at intermediate scales and directed to- 
wards the low- and high-frequency limits, respectively. 
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The respective fluxes were found to be consistent witfi 
the particular scaling exponents of the momentum dis- 
tribution in the respective regimes. Our main result is 
the identification of the nonthermal fixed-point with the 
appearance of topological excitations in the system. We 
could successfully employ statistical models of point vor- 
tices in two dimensions and vortex rings in three dimen- 
sions to interpret the IR scaling exponent of the spectra. 
Moreover, we observed in our simulations, during the fi- 
nal thermalization stage, decay of the turbulence scal- 
ing into a scaling derived from the assumption of vortex- 
antivortex pairing. To characterize further the bimodal 
power spectra developed close to the nonthermal fixed 
point we have analysed the decomposition of the overall 
flow pattern into compressible and incompressible com- 
ponents. We found that the IR power spectra of under- 
lying compressible excitations suggest an understanding 
in terms of acoustic turbulence on top of the vorticity- 
bearing quasicondensate. Finally, to compare with the 
results of recent experiments and analytical predictions, 
the velocity- field probability distribution as well as the 
velocity statistics of individual vortices has been stud- 
ied. This observable is of great interest, as it can be 
used to experimentally distinguish between classical and 
quantum turbulence. 

The connection of QT phenomena with ab-initio dy- 
namical field theoretic methods points a way to unified 
analytical studies of turbulence. Moreover, it provides 
hints of how the proposed nonthermal fixed points in 
relativistic systems [6-8, 10, 14] are realized in nature. 
A manifestation of the approach of a nonthermal fixed- 
point in terms of quasi-solitary pattern formation and 
charge separation in a relativistic scalar model has re- 
cently been reported in Ref. [17]. 

Experimental studies of universal phenomena in 
nonequilibrium dynamics of ultracold atoms have great 
potential since universal effects do not depend signifi- 
cantly on initial conditions and details of the system. 
Following this idea, our numerical protocol was chosen 
such that an experimental verification of our findings is 
within reach of present-day cold-atom experiments. The 
study of turbulence in ultracold gases may have great 
impact on many other fields of physics. Prominent ex- 
amples are strongly correlated nuclear matter produced 
in heavy-ion collisions and early-universe cosmology. 
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Appendix A: Vortical excitations in a superfiuid 

1. Single vortex in d = 2 dimensions 

A singly quantized vortex with circulation k = ±1 
in two dimensions is described, in polar coordinates, 
by the solution (f){r, ip) = f(r) exp{iK(p} of the Gross- 
Pitaevskii equation (10), where /(r) is real and ap- 
proaches the square root of the bulk density ribuik for 
large distances r from the vortex core where /(O) = 0. 
A vortex is a stationary solution of Eq. (10), evolving as 

1 /2 

4>{r,Lp,t) = 4>{r,Lp,Q)eyi■p{-i^it} with = 5'Vuik- The 
velocity field of the vortex is 



(Al) 



The Fourier transform of this field can be conveniently 
evaluated by the help of the first-order Bessel function, 
J d(p cos(v3) exp{iacosiy9} = 27r«Ji(a), as 



m 



k 



(A2) 



2. Straight vortex line in d = 3 

Based on Eq. (31), we calculate the vorticity of two 
straight vortex lines in three dimensions, with circulation 
Ki, i = 1,2, placed at positions (0, —yo) and (0, yo) in the 
x-y plane and relate it to the momentum spectrum, see 
Eq. (37). The lines are parametrized by 



3,;(x) = (0, (-l)Vo,K»T-) 



(A3) 



with r G [— c», oo]. In Fourier space, the vorticity density, 
defined in Eq. (29), reads 



(A4) 



The vorticity spectrum follows from Eq. (31), 

|a;(k)|2 = 2S^{k,) [1 4- KiK2COs(2feyyo)] • (A5) 

Hence, for co- rotating vortex lines (ki = K2) 

\u:{k)\^ = 4S\k,)cos\kyyo) , (A6) 

while for counter-rotating vortex lines (ki = —K2) one 
obtains 



|a;(k)|2 =4<52(fc,)sin2(fc^yo). 



(A7) 
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FIG. 21: Radial momentum distribution of two straight vor- 
tex lines of circulation ki = — K2 = 1 on a A^^ = 1024'^ lat- 
tice. The distance scale kyg is indicated. Note the double- 
logarithmic scale. 

In order to be able to take the angle average over we 
regularize the delta distribution, giving it a finite width 
A = l/Lz- The angle average of Eq. (A6) then reads 

\u}{k)\'^ = J dipkduSA{kuy 

X 4 cos^ {kyo \/l — v?cosLpk ) 
= 87r2i,fc-i[l + Jo(2fcyo)], (AS) 

where Lz is the length of the system in z-direction and 
u — cos{0). One obtains the scaling behaviour 

\u}{k)\^ = WTT^Lzk-^ + 0{k), (A9) 

which is the scaling of a single line. With Eq. (37), we 
conclude n{k) k~^. Following the same reasoning, the 
angle average of Eq. (A7) is given as 

|a;(fc)|2 = 87r2i,fc-i[l- Jo(2fcyo)]. (AlO) 

For kyo < 1 

|a;(A:)|2 = Sn^ L^y^k + Oik^) , (All) 

which gives n{k) ^ k^^. For kyo ^ 1 

|u;(fc)p Sn^Lzk-^ (A12) 

which gives n{k) ^ k~^. These scalings are confirmed in 
Fig. 21. 

3. Ensembles of parallel vortex lines in d = 3 

In Sects. IIIB and III C, we have discussed the scaling 
of the momentum distribution for ensembles of indepen- 
dent and pair correlated vortices in two dimensions, re- 
spectively. In this appendix we generalize the scaling de- 
rived above for straight vortex lines to the case of many 
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FIG. 22: Radial momentum distribution of a vortex ring of 
radius r on a A^^ = 1024^ lattice. The radius scale kr is 
indicated. Note the double-logarithmic scale. 

such lines in d = 3. We derive an expression for the 
momentum spectrum of a system of M straight vortex 
lines oriented along the z direction. The average vortic- 
ity squared is given by Eq. (35) where are vectors in 
the two-dimensional x-y plane pointing to the i-th vortex 
line with circulation Ki — ±1. The vorticity of the single 
line is 

a>i(k) = J d^xe^'^'^y dr e^(5(KjTe^ - x) , 

= K,Sikz)e,. (A13) 

Hence, Eq. (35) has the same form as Eq. (15) for d — 2, 

M 

up to a 5'^{kz) term, which arises from the infinite extent 
of the vortex lines in z-direction. As in App. A 2, we 
regularize the delta distributions, and averaging (A14) 
over solid angles yields 

M 

(l'^(fc)l') - -^Y.^^^^>^=J^{k\■R,^B.,\)), (A15) 

id 

which is analogous to Eq. (16) in two dimensions. As 
the position of the i-th vortex line is determined by Ri, 
statistical averaging of (A15) can be done in the same 
way as in d = 2. 



4. Vortex ring in d = 3 

A vortex ring with radius r lying in the x-y plane can 
be parametrized as 

s(x) = (r cos (/3, r sin (/?, 0) (A16) 
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with (fi G [0, 2tt]. In Fourier space, similar steps as above 
lead to the angle-avcragcd vorticity density 

\Lj{k)\^ = J dOk |Ji(fcrsin0fe)|2 . (A17) 

Evaluating the integral numerically shows that the vor- 
ticity scales like ~ for kr <^ 1 and ^ k^^ for fcr ^ 1, 
corresponding to n{k) ^ and n{k) ^ k~^. This 

is shown in Fig. 22. Both scalings are also found in 
the momentum spectrum of a vortex ellipse, discussed 
in Sect. HID. 



Appendix B: Semiclassical field simulations 

In this paper, a dilute superfluid gas of Bosons of mass 
m is considered, with contact interactions quantified by 
a coupling constant g. Its dynamics is described, in 
the classical wave limit, by the Gross-Pitaevskii equa- 
tion (GPE) (10). We use units where H = 1. We 
consider the gas to be contained in a box of size L'^, 
d = 2, 3, with periodic boundary conditions. In d = 3, 
the coupling constant g is related to the s-wave scat- 
tering length a hy g = gsu = 4Tra/m. In c? = 2, 
one has g = 520 = ^(47r/rn)[ln(/ima|j3/4)]~^ where 
fi is the chemical potential and 020 is a scattering 
length in two dimensions which, for a gas of hard-spheres 
of radius a is given by a2u = ae^, with the Euler- 
Mascheroni constant 7 ~ 0.577 [78-80]. For a two- 
dimensional gas created by trapping a three-dimensional 
one tightly in one dimension, with harmonic-oscillator 
length Iz, the effective 2D scattering length is given by 
a2D = 4/^(7r/B)i/2exp{-V^ZJa3D}, where B ~ 0.915 
[81]. Dilutencss implies that a <^ I, the interparticlc 
spacing I = n~^^'^ being determined by the mean density 
n = N/L"^. 

The initial values for the real and imaginary parts 
of the field (/)(k, 0) are randomly chosen from a Gaus- 
sian distribution with width 1/2, centered around 
\/n(k,0) exp{i(^(k,0)}, where n{k,t) = (0t(k, i)0(k, t)) 
is the occupation number at time t and (/'(k, 0) is a ran- 
dom phase angle. Correlation functions including n(k, t) 
are obtained by averaging over many trajectories. To in- 
duce transport from small to large wave numbers, only a 
few modes near k = are chosen to be macroscopically 
occupied at the initial time, ri(k, 0) 3> 1. Such an ini- 
tial state can be prepared, e.g., by Bragg scattering of 
photons from a Bose-Einstein condensate. 

Our numerical simulations are performed on space- 
time lattices with side lengths L — NgUg, with lattice 
spacing and in total Nf gridpoints. We use periodic 
boundary conditions and Ng G {256, 512} for d = 2 and 
Ns e {128,256,512} for d = 3. Eq. (10) is written in 
terms of the dimensionless variables 'g = 2mgal~^^, t = 
t/(2mag) and V'nCO = V'n-\/a?Gxp(2it). All quantities in 
the figures are either dimensionless or shown in lattice 
units, with the length unit given by a^. The dimension- 
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FIG. 23: (Color online) Normalized occupation number differ- 
ence \n{k)—ns{k)\/n{k) as a function of the radial momentum 
k at different times, for d = 2 (left panel) and d = 3 (right 
panel). Parameters are the same as in Fig. 10, for the run 
in d = 2 dimensions and in Fig. 11, for d = 3. Note the 
logarithmic scale. 

less lattice momenta are k = [^^^j^ 4sin^(fci/2)]^/^, k = 
2™/iV„ n = (ni,...,nd), n, = -NJ2, N,/2. 

In order to relate our simulations to a typical situa- 
tion in experiment, we give parameters for Rb-87. We 
estimate the total energy of the gas in equilibrium to 
be given by the interaction energy which, by equipar- 
tition is related to the temperature, 2Eint = L'^gn^ ~ 
Nfafgn^ = N^ksT. At the transition to degeneracy the 
thermal de Broglie wave length AdB = y/2TTfi? /rnksT 
is of the order of the interparticlc spacing and hence 
ksTdcg ^ 2T:h^ii?/'^ /m. Inserting this into the expres- 
sion for the interaction energy allows to express the lat- 
tice spacing in terms of the density and the scattering 
length as 



«s = 7;-V^lii4 - ln(/Ltma2D) (d = 2), 
Zn 

a, = (n'^/3^)-i/3 (d = 3). (Bl) 

Inserting parameters for a typical Rb-87 experiment in 
d = 3, viz., a ~ 5 nm, m ~ 1.4 x 10~^^kg, n = 
10^" m^'^ we obtain a,, = 1 ^m. Hence, time scale as 
t/t = 2mal/h'^ = 3 X 10^^ s. Our lattice time typically 
runs until t ~ 1000, which therefore represents a realistic 
observation timescale in experiments. Typical parame- 
ter choices in our simulations are N = 10^, = 128, 
g = 4 X 10"^ for d = 3 and iV = 4 X 10^ A^,, = 512, 
g = 3 X 10"^ for d = 2. 



Appendix C: Decomposition of w 

In Sect. VA, we have defined occupation numbers 
corresponding to the incompressible, compressible and 
quantum pressure components of the kinetic energy. Nu- 
merically, the resulting spectra ni(fc), nc(fc), and nq(fc) 
add up to ns{k) which is equal to the full single-particle 
spectrum n(fc) up to small corrections. In particular, the 
scaling behavior emerges as the same for both, ns{k) and 
ri(fc). To make this more explicit, we calculate k'^n{k) in 
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terms of generalized velocities wg, S S {v,q}, 

k^n{k) J-(wve*'^)(fc) J-(wve''^)* (fc) ) 

+ 23fJ( iJ'(wve"^)(/s) J'(wqe*'^)*(fc) ) 

+ ( J-(wqe*'^)(fc)^(wqe^^)*(fc)) . (CI) 

3? denotes the real part and J-' the Fourier transform. 
The middle term vanishes since the expectation value is 
purely imaginary due to isotropy. It follows that 

k'^n{k) 

+ ( * ^e'nm (wq * He'nr(k) ) . (C2) 



Since phase and phase velocity are expected not to be 
correlated, the main contribution to the 4-point correla- 
tion functions is assumed to arise from terms of the form 
(W5(|p-k|)w5(|q-k|))(j-(exp{*^})(p)^(exp{*^})*(g)), 
6 G {v,q}. In the superfluid regime, the phase of the field 
is a slowly varying function in position space. Therefore, 
J^(exp{i(^})*((jf)-terms are strongly peaked at zero mo- 
mentum, effectively acting as regularized delta functions 
under the convolution. Hence, 

fcVfc)=^(|wv(fc)|'> + (K(fc)|2). (C3) 

This approximative identity is supported by Fig. 23. 



[1] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov 

(Cambridge University Press, Cambridge, UK, 1995). [24 
[2] L. F. Richardson, Proc. Roy. Soc. Lend. A 97, 354 (1920). [25 
[3] A. N. Kolmogorov, Dokl. Akad. Nauk. USSR 30, 299 

(1941)[Sov. Phys. Dokl. 10, 734 (1968)]; reprinted in 

Proc. Roy Soc. bond. A, 434, 9 (1991). [26 
[4] A. M. Obukhov, Izv. Akad. Nauk S.S.S.R., Ser. Geogr. 

Geofiz. 5, 453 (1941). [27 
[5] R. Micha and I. I. Tkachev, Phys. Rev. Lett. 90, 121301 

(2003). [28 
[6] J. Berges, A. Rothkopf, and J. Schmidt, Phys. Rev. Lett. 

101, 041603 (2008). [29 
[7] J. Berges and G. Hoffmeister, Nucl. Phys. B813, 383 

(2009). [30 
[8] J. Berges, S. Scheffler, and D. Sexty, Phys. Lett. B681, 

362 (2009). [31 
[9] C. Scheppach, J. Berges, and T. Gasenzer, Phys. Rev. A 

81, 033611 (2010). [32 
10] J. Berges and D. Sexty, Phys. Rev. D 83, 085004 (2011). 
11] P. B. Arnold and G. D. Moore, Phys. Rev. D 73, 025006 [33 

(2006). 

2] P. B. Arnold and G. D. Moore, Phys. Rev. D 73, 025013 [34 
(2006). [35 

3] A. H. Mueller, A. L Shoshi, and S. M. H. Wong, Nucl. 

Phys. B760, 145 (2007); [36 

4] M. Carrington and A. Rebhan, Eur. Phys. J. C71, 1787 

(2011) . [37 
5] K. Fukushima and F. Gelis, Nucl. Phys. A 874, 108 

(2012) . 

6] B. Nowak, D. Sexty, and T. Gasenzer, Phys. Rev. B 84, [38 

020506(R) (2011). 
7] T. Gasenzer, B. Nowak, and D. Sexty, Phys. Lett. B 710, [39 
500 (2012). 

8] W. P. Halperin and M. Tsubota, eds.. Progress in Low [40 
Temperature Physics Vol. 16: Quantum Turbulence (El- 
sevier, Amsterdam, 2008). [41 
9] R. J. Donnelly, Quantized Vortices in Liquid He II (CUP, 

Cambridge, 1991). 
0] T.-L. Horng, C.-H. Hsueh, and S.-C. Gou, Phys. Rev. A [42 

77, 063625 (2008). 
1] T.-L. Horng, C.-H. Hsueh, S.-W. Su, Y.-M. Kao, and [43 

S.-C. Gou, Phys. Rev. A 80, 023618 (2009). 
2] V. Yukalov, Las. Phys. Lett. 7, 467 (2010). [44 
3] C. J. Foster, P. B. Blakie, and M. J. Davis, Phys. Rev. 



A 81, 023623 (2010). 

L. Onsager, Nuovo Cim. Suppl. 6, 279 (1949). 
R. P. Feynman, in Progress in Low Temperature Physics 
Vol 1., edited by C. J. Gorter (North Holland, Amster- 
dam, 1955), p. 17. 

J. Maurer, and P. Tabeling, Europhys. Lett. 43, 29 
(1998). 

S. R. Stalp, L. Skrbek, and R. J. Donnelly, Phys. Rev. 
Lett. 82, 4831 (1999). 

W. F. Vinen and J. J. Niemela, J. Low Temp. Phys. 128, 
167 (2002). 

C. Nore, M. Abid, and M. E. Bracket, Phys. Rev. Lett. 
78, 3896 (1997). 

C. Nore, M. Abid, and M. E. Brachet, Phys. Fluids 9, 
2644 (1997). 

T. Araki, M. Tsubota, and S. K. Nemirovskii, Phys. Rev. 
Lett. 89, 145301 (2002). 

M. Kobayashi and M. Tsubota, Phys. Rev. Lett. 94, 
065302 (2005). 

M. Kobayashi and M. Tsubota, J. Phys. Soc. Jpn. 74, 
3248 (2005). 

M. Tsubota, J. Phys. Soc. Jpn. 77, 111006 (2008). 

G. Krstulovic and M. Brachet, Phys. Rev. E 83, 066311 

(2011). 

E. Levich and V. Yakhot, J. Phys. A: Math. Gen. 11, 
2237 (1978). 

Y. Kagan, B. V. Svistunov, and G. V. Shiyapnikov, Zh. 
Eksp. Teor. Fiz. 101, 528 (1992) [Sov. Phys. JETP 74, 
279 (1992)]. 

Y. Kagan and B. V. Svistunov, Zh. Eksp. Teor. Fiz. 105, 

353 (1994) [Sov. Phys. JETP 78, 187 (1994)]. 

Y. Kagan and B. V. Svistunov, Phys. Rev. Lett. 79, 3331 

(1997). 

N. G. Berloff and B. V. Svistunov, Phys. Rev. A 66, 
013603 (2002). 

B. Svistunov, in Quantized Vortex Dynamics and Super- 
fluid Turbulence, edited by C. Barenghi, R. Donnelly, and 
W. Vinen (Springer, Berlin, 2001). 

E. V. Kozik and B. V. Svistunov, J. Low Temp. Phys. 
156, 215 (2009). 

C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, 
M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008). 
E. A. L. Henn, J. A. Seman, G. Roati, K. M. F. Mag- 
alhaes, and V. S. Bagnato, Phys. Rev. Lett. 103, 045301 



21 



(2009). 

J. A. Seman, E. A. L. Henn, R. F. Shiozaki, G. Roati, [64] 

F. J. Poveda-Cuevas, K. M. F. Magalhaes, V. I. Yukalov, 

M. Tsubota, M. Kobayashi, K. Kasamatsu, and V. S. [65 
Bagnato, Laser Physics Letters 8, 691 (2011). [66 
V. E. Zakharov, V. S. LVov, and G. Falkovich, Kol- [67 
mogorov Spectra of Turbulence I: Wave Turbulence 
( Springer- Verlag, Berlin, 1992). [68 
S. Nazarenko, Wave turbulence, no. 825 in Lecture Notes 
in Physics (Springer, Heidelberg, 2011). [69 
A. Newell, Ann. Rev. Fluid Mech. 43 (2011). 
M. S. Paoletti, M. E. Fisher, K. R. Sreenivasan, and D. P. [70 
Lathrop, Phys. Rev. Lett. 101, 154501 (2008). 
J. M. Luttinger and J. C. Ward, Phys. Rev. 118, 1417 [71 
(1960). 

G. Baym, Phys. Rev. 127, 1391 (1962). [72 
J. M. Cornwall, R. Jackiw, and E. Tomboulis, Phys. Rev. 
D 10, 2428 (1974). [73 
J. Berges, Nucl. Phys. A699, 847 (2002). [74 
For videos of the evolution see: http://www.thphys.uni- [75 
heidelberg.de / ~smp/gasenzer / videos/boseqt .html 
G. Aarts, D. Ahrensmeier, R. Baier, J. Berges, and J. Ser- [76 
reau, Phys. Rev. D 66, 045008 (2002). 

G. Eyink and K. Sreenivasan, Rev. Mod. Phys. 78, 87 
(2006). [77 
P. Chavanis, Statistical mechanics of two-dimensional 
vortices and stellar systems (Springer, Berlin, 2002). [78 
V. L. Berdichevsky, Phys. Rev. E 51, 4432 (1995). [79 
A. C. White, C. F. Barenghi, N. P. Proukakis, A. J. Youd, 
and D. H. Wacks, Phys. Rev. Lett. 104, 075301 (2010). [80 
E. A. Novikov, Zh. Eksp. Teor. Fiz. 68, 1868 (1975); [Sov. 
Phys. JETP 41, 937 (1975)]. [81 
C. Lin, Proc. Nat. Acad. Sci. 27, 570 (1941). 
S. K. Nemirovskii, Phys. Rev. B 57, 5972 (1998). 
S. K. Nemirovskii, M. Tsubota, and T. Araki, J. Low 



Temp. Phys. 126, 1535 (2002). 

A. Chorin and J. Akao, Physica D: Nonlin. Phen. 52, 403 
(1991). 

V. L. Berdichevsky, Phys. Rev. E 57, 2885 (1998). 
V. L. Berdichevsky, Int. J. Eng. Sci. 40, 123 (2002). 
S. Nemirovskii and L. Kondaurova, J. Low Temp. Phys. 
156, 182 (2009). 

J. B. Weiss, A. Provenzale, and J. C. McWilliams, Phys. 
Fluids 10, 1929 (1998). 

I. A. Min, I. Mezic, and A. Leonard, Phys. Fluids 8, 1169 
(1996). 

P.-H. Chavanis and C. Sire, Physics of Fluids 13, 1904 
(2001). 

S. Nazarenko and M. Onorato, Physica D: Nonlin. Phen. 
219, 1 (2006). 

R. Numasato, M. Tsubota, and V. S. L'vov, Phys. Rev. 
A 81, 063630 (2010). 

S. Khlebnikov, Phys. Rev. A 66, 063606 (2002). 
K. E. Daniels and E. Bodenschatz, Chaos 13, 55 (2003). 
D. Proment, S. Nazarenko, and M. Onorato, Physica D: 
Nonlin. Phen. 241, 304 (2012). 

M. Leadbeater, T. Winiecki, D. C. Samuels, 

C. F. Barenghi, and C. S. Adams, Phys. Rev. Lett. 86, 
1410 (2001). 

N. G. Berloff and A. J. Youd, Phys. Rev. Lett. 99, 145301 
(2007). 

M. Schick, Phys. Rev. A 3, 1067 (1971). 

D. S. Fisher and P. C. Hohenberg, Phys. Rev. B 37, 4936 
(1988). 

M. D. Lee, S. A. Morgan, M. J. Davis, and K. Burnett, 
Phys. Rev. A 65, 043617 (2002). 

D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, 
Phys. Rev. Lett. 84, 2551 (2000). 



